#region License
// Copyright 2007-2008 Rory Plaire (codekaizen@gmail.com)
/*
 * This license governs use of the accompanying software. If you use the software, you
 * accept this license. If you do not accept the license, do not use the software.
 * 
 * 1. Definitions
 * The terms "reproduce," "reproduction," "derivative works," and "distribution" have the
 * same meaning here as under U.S. copyright law.
 * A "contribution" is the original software, or any additions or changes to the software.
 * A "contributor" is any person that distributes its contribution under this license.
 * "Licensed patents" are a contributor's patent claims that read directly on its contribution.
 * 
 * 2. Grant of Rights
 * (A) Copyright Grant- Subject to the terms of this license, including the license conditions 
 *     and limitations in section 3, each contributor grants you a non-exclusive, worldwide, 
 *     royalty-free copyright license to reproduce its contribution, prepare derivative works 
 *     of its contribution, and distribute its contribution or any derivative works that you create.
 * (B) Patent Grant- Subject to the terms of this license, including the license conditions and 
 *     limitations in section 3, each contributor grants you a non-exclusive, worldwide, royalty-free 
 *     license under its licensed patents to make, have made, use, sell, offer for sale, import, 
 *     and/or otherwise dispose of its contribution in the software or derivative works of the 
 *     contribution in the software.
 * 
 * 3. Conditions and Limitations
 * (A) No Trademark License- This license does not grant you rights to use any contributors' 
 *     name, logo, or trademarks.
 * (B) If you bring a patent claim against any contributor over patents that you claim are 
 *     infringed by the software, your patent license from such contributor to the software 
 *     ends automatically.
 * (C) If you distribute any portion of the software, you must retain all copyright, patent, 
 *     trademark, and attribution notices that are present in the software.
 * (D) If you distribute any portion of the software in source code form, you may do so only under 
 *     this license by including a complete copy of this license with your distribution. If you 
 *     distribute any portion of the software in compiled or object code form, you may only do so 
 *     under a license that complies with this license.
 * (E) The software is licensed "as-is." You bear the risk of using it. The contributors give no 
 *     express warranties, guarantees or conditions. You may have additional consumer rights under 
 *     your local laws which this license cannot change. To the extent permitted under your local laws, 
 *     the contributors exclude the implied warranties of merchantability, fitness for a particular 
 *     purpose and non-infringement.
 * 
 */
#endregion

// Based initially on JAMA : A Java Matrix Package by 
// National Institute of Standards and Technology (NIST) 
// and The MathWorks, a public-domain work.

using System;
using NPack.Numeric;
using NPack.Numeric.Interfaces;

namespace NPack
{
    /// <summary>
    /// LU decomposition of a rectangular matrix.
    /// </summary>
    /// <typeparam name="T">The type of the matrix element.</typeparam>
    /// <remarks>
    /// <para>
    /// For an m-by-n Matrix{T} <c>A</c> with m &gt;= n, the LU decomposition is an m-by-n
    /// unit lower triangular Matrix{T} <c>L</c>, an n-by-n upper triangular Matrix{T} <c>U</c>,
    /// and a permutation vector <c>Pivot</c> of length m so that <c>Pivot * A = L * U</c>.
    /// If m &lt; n, then <c>L</c> is m-by-m and <c>U</c> is m-by-n.
    /// </para>
    /// <para>
    /// The LU decompostion with pivoting always exists, even if the Matrix{T} is
    /// singular, so the constructor will never fail.  The primary use of the
    /// LU decomposition is in the solution of square systems of simultaneous
    /// linear equations. There is no solution, or infinite solutions if 
    /// <see cref="IsNonsingular"/> returns <see langword="false"/>.
    /// </para>
    /// </remarks>
    public class LUDecomposition<TComponent, TMatrix> : ILUDecomposition<TComponent, TMatrix>
        //where TComponent : IEquatable<TComponent>, IComparable<TComponent>, IComputable<TComponent>, IConvertible, IFormattable
        where TMatrix : IMatrix<TComponent, TMatrix>, IEquatable<TMatrix>, IComparable<TMatrix>, IComputable<TMatrix>
    {
        private readonly TMatrix _inputMatrix;
        private readonly TMatrix _lu;
        private readonly Int32 _pivotSign;
        private readonly Int32[] _rowPivotVector;
        private readonly IMatrixFactory<TComponent, TMatrix> _matrixFactory;

        /// <summary>
        /// Create a new LU decomposition.
        /// </summary>	
        public LUDecomposition(TMatrix value, IMatrixFactory<TComponent, TMatrix> matrixFactory)
        {
            if (matrixFactory == null)
            {
                throw new ArgumentNullException("matrixFactory");
            }

            if (Equals(value, default(TMatrix)))
            {
                throw new ArgumentNullException("value");
            }

            _matrixFactory = matrixFactory;

            _inputMatrix = (TMatrix)value.Clone();
            _lu = (TMatrix)value.Clone();

            Int32 rowCount = _lu.RowCount;
            Int32 columnCount = _lu.ColumnCount;
            _rowPivotVector = new Int32[rowCount];

            for (Int32 i = 0; i < rowCount; i++)
            {
                _rowPivotVector[i] = i;
            }

            _pivotSign = 1;

            // Outer loop.
            for (Int32 j = 0; j < columnCount; j++)
            {
                // Apply previous transformations.
                for (Int32 i = 0; i < rowCount; i++)
                {
                    // Most of the time is spent in the following dot product.
                    Int32 kmax = Math.Min(i, j);
                    TComponent dotProduct = default(TComponent);

                    for (Int32 k = 0; k < kmax; k++)
                    {
                        dotProduct = dotProduct + ((Real<TComponent>)_lu[i, k] * _lu[k, j]);
                    }

                    _lu[i, j] = (Real<TComponent>)_lu[i, j] - dotProduct;
                }

                // Find pivot and exchange if necessary.
                Int32 pivotIndex = j;

                for (Int32 i = j + 1; i < rowCount; i++)
                {
                    if (((Real<TComponent>)_lu[i, j]).Abs() > ((Real<TComponent>)_lu[i, pivotIndex]).Abs())
                    {
                        pivotIndex = i;
                    }
                }

                if (pivotIndex != j)
                {
                    for (Int32 k = 0; k < columnCount; k++)
                    {
                        TComponent t = _lu[pivotIndex, k];
                        _lu[pivotIndex, k] = _lu[j, k];
                        _lu[j, k] = t;
                    }

                    Int32 v = _rowPivotVector[pivotIndex];
                    _rowPivotVector[pivotIndex] = _rowPivotVector[j];
                    _rowPivotVector[j] = v;

                    _pivotSign = -_pivotSign;
                }

                // Compute multipliers.

                if (j < rowCount && !_lu[j, j].Equals(default(TComponent)))
                {
                    for (Int32 i = j + 1; i < rowCount; i++)
                    {
                        _lu[i, j] = (Real<TComponent>)_lu[i, j] / _lu[j, j];
                    }
                }
            }
        }

        #region IFactorization<T> Members

        public TMatrix InputMatrix
        {
            get { return (TMatrix)_inputMatrix.Clone(); }
        }

        #endregion

        /// <summary>
        /// Gets true if the <see cref="Matrix{T}"/> is non-singular, 
        /// false if singular (non-invertible).
        /// </summary>
        public Boolean IsNonsingular
        {
            get
            {
                for (Int32 j = 0; j < _lu.ColumnCount; j++)
                {
                    if (_lu[j, j].Equals(default(TComponent)))
                    {
                        return false;
                    }
                }

                return true;
            }
        }

        /// <summary>
        /// Gets the determinant of the matrix.
        /// </summary>
        /// <exception cref="InvalidOperationException">Thrown if the matrix used to compute the factorization is non-square.</exception>
        public Double Determinant
        {
            get
            {
                if (_lu.RowCount != _lu.ColumnCount)
                {
                    throw new InvalidOperationException("Matrix<T> must be square.");
                }

                Double determinant = _pivotSign;

                for (Int32 j = 0; j < _lu.ColumnCount; j++)
                {
                    determinant *= Convert.ToDouble(_lu[j, j]);
                }
                return determinant;
            }
        }

        /// <summary>
        /// Gets the lower triangular matrix <c>L</c> such that 
        /// <c>A * <see cref="Pivot"/> = L * <see cref="U"/></c>.
        /// </summary>
        public TMatrix L
        {
            get
            {
                Int32 rows = _lu.RowCount;
                Int32 columns = _lu.ColumnCount;
                TMatrix l = _matrixFactory.CreateMatrix(_lu.Format, rows, columns);

                TComponent one = default(Real<TComponent>).One;

                for (Int32 i = 0; i < rows; i++)
                {
                    for (Int32 j = 0; j < columns; j++)
                    {
                        if (i > j)
                        {
                            l[i, j] = _lu[i, j];
                        }
                        else if (i == j)
                        {
                            l[i, j] = one;
                        }
                        else
                        {
                            l[i, j] = default(TComponent);
                        }
                    }
                }

                return l;
            }
        }

        /// <summary>
        /// Gets the lower triangular matrix <c>L</c> such that 
        /// <c>A * <see cref="Pivot"/> = <see cref="L"/> * U</c>.
        /// </summary>
        public TMatrix U
        {
            get
            {
                Int32 rows = _lu.RowCount;
                Int32 columns = _lu.ColumnCount;
                TMatrix u = _matrixFactory.CreateMatrix(_lu.Format, rows, columns);

                for (Int32 i = 0; i < rows; i++)
                {
                    for (Int32 j = 0; j < columns; j++)
                    {
                        if (i <= j)
                        {
                            u[i, j] = _lu[i, j];
                        }
                        else
                        {
                            u[i, j] = default(TComponent);
                        }
                    }
                }

                return u;
            }
        }

        /// <summary>
        /// Gets the pivot permutation matrix.
        /// </summary>
        public TMatrix Pivot
        {
            get
            {
                TMatrix pivot = ((IHasZero<TMatrix>)_lu).Zero;

                for (Int32 i = 0; i < _lu.RowCount; i++)
                {
                    pivot[i, _rowPivotVector[i]] = default(Real<TComponent>).One;
                }

                return pivot;
            }
        }

        /// <summary>
        /// Solves a set of equation systems of type <c>A * X = B</c>.
        /// </summary>
        /// <param name="value">
        /// Right hand side <see cref="Matrix{T}"/> with as many rows as 
        /// <c>A</c> and any number of columns.
        /// </param>
        /// <returns>
        /// <see cref="Matrix{T}"/> <c>X</c> so that <c>L * U * X = B</c>.
        /// </returns>
        public TMatrix Solve(TMatrix value)
        {
            if (Equals(value, default(TMatrix)))
            {
                throw new ArgumentNullException("value");
            }

            if (value.RowCount != _lu.RowCount)
            {
                throw new ArgumentException("Invalid Matrix<T> dimensions.", "value");
            }

            if (!IsNonsingular)
            {
                throw new InvalidOperationException("Matrix<T> is singular");
            }

            // Copy right hand side with pivoting
            Int32 count = value.ColumnCount;
            TMatrix pivot = (TMatrix)value.GetMatrix(_rowPivotVector, 0, count - 1);

            Int32 columnCount = _lu.ColumnCount;

            // Solve L*Y = B(pivot,:)
            for (Int32 k = 0; k < columnCount; k++)
            {
                for (Int32 i = k + 1; i < columnCount; i++)
                {
                    for (Int32 j = 0; j < count; j++)
                    {
                        pivot[i, j] = pivot[i, j] - ((Real<TComponent>)pivot[k, j] * _lu[i, k]);
                    }
                }
            }

            // Solve U*X = Y;
            for (Int32 k = columnCount - 1; k >= 0; k--)
            {
                for (Int32 j = 0; j < count; j++)
                {
                    pivot[k, j] = (Real<TComponent>)pivot[k, j] / _lu[k, k];
                }

                for (Int32 i = 0; i < k; i++)
                {
                    for (Int32 j = 0; j < count; j++)
                    {
                        pivot[i, j] = pivot[i, j] - ((Real<TComponent>)pivot[k, j] * _lu[i, k]);
                    }
                }
            }

            return pivot;
        }
    }
}